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Abstract 

We present numerical evidence for the fact that the damage spreading transition in the 
Domany-Kinzel automaton found by Martins et al. is in the same universality class as 
directed percolation. We conjecture that also other damage spreading transitions should 
be in this universality class, unless they coincide with other transitions (as in the Ising 
model with Glauber dynamics) and provided the probability for a locally damaged state 
to become healed is not zero. 



PACS numbers: 05.20, 05.40, 05.50, 36.20.E, 61.40, 64.70 
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Among all critical phenomena, directed percolation (DP) is maybe that which has 
been associated with the most wide variety of phenomena. 

First there are interpretations where the preferred direction is a spatial direction. 
This was of course proposed to apply to material and charge transport in disor- 
dered media under the influences of external forces. Also, it should model the 
propagation of epidemics and forest fires under some directional bias, e.g. strong 
wind. 

More interesting are interpretations where the preferred direction is time. Here, 
the primary interpretation is as an epidemic without immunization, the so called 
"contact process" p] or the "simple epidemic" ||. 

But these are by no means all possible applications. A very early application (even 
if it took rather long until it was understood as such || [|]) was to "reggeon field 
theory", a theory for ultrarelativistic particle collisions popular in the 70's 01. 
Here, the preferred direction is that of "rapidity" , while the directions transverse 
to it are provided by the impact parameter plane. This connection is interest- 
ing since it was through it that first precise estimates of critical exponents and 
amplitudes were obtained for DP [|5| . 

Another realization of the DP transition occurs in simple models of heterogeneous 
catalysis. The first such model was proposed by Ziff et al. |J (ZGB). The simula- 
tions of these and subsequent authors indicated that this model was in a different 
universality class, and it is only after some controversies that it is now gener- 
ally accepted to be in the DP universality class J7J. Similar models are invented 
again and again B |9| . Repeatedly they are claimed to be in different universality 



classes, and repeatedly these claims are refuted [10, [II], |12 



In ||13| |14| it was proposed that the universality class of DP contains all continuous 
transitions from a "dead" or "absorbing" state to an "active" one with a single 
scalar order parameter, provided the dead state is not degenerate (and provided 
some technical points are fulfilled: short range interactions both in space and 
time, nonvanishing probability for any active state to die locally, translational 
invariance [absence of 'frozen' randomness], and absence of multicritical points). 
It seems fair to say that there is now ample evidence for this proposal. It predicts, 
e.g. immediately that the ZGB model is in this universality class. A rather subtle 
question is whether also chaotic systems where the random noise is replaced by 



deterministic chaos are in the same class [15, 16 



As far as I am aware of, no model with non-fluctuating absorbing state and a multi- 
component order parameter has ever been studied in the literature. Notice that 
this has to be distinguished from models for which some mean field approximation 
has a mult i- component order parameter. Such models are quite common (e.g. the 
Bethe-Peierls approximation of the Ising model, or the mean field approximation 
of ZGB), and it was just the study of such a model which had led to the conjecture 
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in [14]. A supposed generalization [17] of the above conjecture is thus already fully 
contained in the original conjecture of []I1| . 

A more interesting question is what happens if the dead state is degenerate. 
Counter examples with twofold degeneracy were studied in |1| [19], p0[ . They 
involve conservation laws which prevent some active states from dying, making it 
thus immediately clear that any transition — if it occurs at all — has to be in a 
different universality class. 

But the main open problem is whether models can be generically in the DP 
class if they have an absorbing state with positive entropy (for obvious reason, we 
prefer not to call it "dead" in this case). One might conjecture that such a state is 
essentially unique on a coarse scale, provided its evolution is ergodic and mixing - 
and provided it does not involve long range correlations (long correlations should 
be entirely due to patches of "active" states). Since only coarse-grained properties 
should influence critical behavior, this would suggest that such transitions are in 



the DP class. This seems contradicted by simulations of some catalysis models |2T 
p2| . But as we pointed out already, systematic errors are often underestimated 
and recent simulations of these and similar models support the conjecture |24 



25], (as explained below, we believe that violations of universality observed in 
some of the latter papers for 'dynamic' properties should be disregarded). 

In this note we propose that there is a rather large and well studied class of tran- 
sitions which are exactly of the latter type, and which are thus all in the DP 
class. These are so called "damaging" transitions. In these models one considers 
two replicas of a stochastic spin system, and lets them evolve with identical re- 
alizations of the stochastic noise. The initial conditions can be either completely 
independent, or one can start with two states which are identical except for a sin- 
gle spin. This single flip is considered as a "damage" , and the question is whether 
this damage will finally heal, so that both replicas converge towards identical 
states — or whether it will spread. If the two states are uncorrelated initially, 
the transition is between a situation where their rescaled Hamming distance (= 
density of damaged sites) stays finite and one where this distance goes to zero. 

More precisely, we propose that such damaging transitions are in the DP class 
if they do not coincide with another transition (since then there would be long 
range correlations in the absorbing state), and if there is no frozen randomness. 
The former applies to the 2-d Ising model with Glauber dynamics, since there the 
damaging transition coincides with the ordinary critical point |27| (the situation 



is less clear in 3 dimensions |]2"q, |29"f). Frozen randomness is involved in damaging 



in spin glasses [pO|, pl| , p^j and in the extensive studies of damaging in Kauffman 
models |33|. |34|. [35| . This should be in the same universality class as DP with frozen 
randomness | 36| , but for the Kauffman models there is a further complication: 
there a damage typically does not heal completely, whence the 'dead' state is not 
absorbing in our sense. We might mention that it was already pointed out that 



damaging in the annealed Kaufmann model is in the DP class ||34||, but this is 
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much more trivial than our present claim. The annealed model can be mapped 
exactly onto DP, which is not the case in general. 



We support our claim with simulations of damaging in the Domany-Kinzel cellular 
automaton (CA) [[37]]. This is a CA with one space and one time dimension, and 
with two states per site: Si = 0,1. Dynamics is defined by the following rule 
involving two real parameters p\ and p 2 (we make a trivial modification which 
slightly simplifies the simulation): 

(i) if Si = and s i+1 = then s' { = 

(ii) if Si XOR Si + i = 1 then s^ — l with probability p\ and = with probability 
1 - pi 

(iii) if Si AND s i+ i = 1 then s[ = 1 with probability p 2 and s[ = with probability 
1 -P2- 

For pi < 1/2 and any p 2 < 1, it is obvious that any state will converge towards 
the dead state . . . 000 .... Actually, this state is an attractor for all values of p\ 
below a critical curve pl(p 2 )- This curve is indicated as curve C in fig.l. To the 
right of C, one has an active state (the dead state still is stationary, but it no 
longer attracts all initial states) with p = (sj) > 0. 

The above conjecture suggests that the transition all along C is in the DP class, 
except at its upper limit point (p±,p 2 ) = (1/2, 1) where the model is a discrete 
time variant of the exactly solvable voter model [0 ("compact directed percola- 
tion" This is supported by all numerical evidence |39| ||| (except for a 



renormalization group analysis and Monte Carlo simulations presented in ^T|; in 
high precision Monte Carlo simulations [|12|] we could not confirm these claims). 
In particular, bond and site DP correspond to p 2 = (2 — p\)p± and p\ = p 2 , 
respectively. 



It was found recently in that the active phase can be further subdivided 



into a phase in which damage does not spread ("healing active phase") and one 
where it does ("chaotic"). The transition between these two phases is indicated 
by curve V in fig.l. It corresponds to pi = pf{p 2 ) where p\ > pi for all p 2 > 0, 
while pf(0) = pl(0) H|. Indeed, one can consider two different variants of the 
damaging process: in the first one uses different random numbers when applying 
rules (ii) and (iii) above, in the second one uses the same. Curve V is computed 
with the second variant. The first variant would give a different curve slightly to 
the left of V Q. 



As pointed out in P0"|, one can describe a pair of replicas by an extended phase 



space with 4 states per site: (00), (01), (10) and (11). Damage spreading corre- 
sponds then to the (directed) percolation of states (10) and (01), while any state 
with (00) and (11) only is healed. Since pf > pi for all p 2 > 0, the healing state has 
positive entropy at the damage spreading transition, and it does not immediately 

1 The very existence of these two variants shows that it is misleading to speak of different 
phases in the Domany-Kinzel CA, as done in [Q. Instead these are different phases for very 
specific algorithms for simulating pairs of such automata. 
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follow from the conjecture of [13, 14] that this transition is in the DP universality 
class. 



To check our conjecture that it is in this class nevertheless, we performed extensive 
simulations at p\ — 1 where both variants coincide. Less extensive runs were made 
at several other values of pi, where we studied both variants. 

We worked on lattices of length L with periodic boundary conditions. To speed up 
our simulations, we simulated 64 lattices simultaneously (we worked on machines 
with 64 bit words) by assigning the k-th bit of the i-th word in an integer array 
of length L to the spin sj in the k-th lattice. The dynamics is then easily 
implemented by standard bit operations. 

To measure the degree of damage in simulations which start with independent 
random initial configurations (thus with half of the sites damaged initially), we 
count the number of set bits in each word. If this number is rii for the i-th word, 
then the number of pairs of lattices which are damaged at site i is (64 — ni)rii. 
The sum of Hamming distances between all 64 x 63/2 = 2016 pairs of lattices is 
thus 

L 

d = J2 (64 - n^m . (1) 

i=i 

For simulations with initial single site damage, this is not possible since we cannot 
build an initial state in which each pair is damaged at only one site. Instead, 
we introduced single site damage only between successive bits, i.e. we initially 
damaged the (2k + l)-st bit (k — 0, 1, 2, . . . 31) in 32 different words, and counted 
how often the (2k + l)-st bit differed from the (2fc)-th one. 

Results from runs with random and independent initial states on lattices of size 
L = 2 22 are presented in fig. 2. There we show the total damage as function of 
time, for p± = 1 and several values of p 2 . At the critical point we expect an 
algebraic decay, corresponding to a straight line in fig. 2. If the transition is in the 
DP class, this decay is governed by an exponent 6 = 0.1596 ± 0.0001 H, HR PJ. 



We see indeed a nearly perfect straight line for p 2 ~ 0.3122. Together with the 
data described below this gives our estimate 

p d 2 = 0.31215 ± 0.00004 for p x = 1 , (2) 

and the exponent extracted from it (0.157± 0.002) is in good agreement with DP. 
Similar results (although somewhat less precise) were obtained for both variants 
of the damage spreading at several values of p x . For pi = 0.85, e.g., they gave 
p d 2 = 0.1957 ± 0.0002 (variant 1) resp. p d 2 = 0.1400 ± 0.0002 (variant 2). In all 
cases 5 was compatible with the DP value. 

It is well known from studies of DP that the exponent (3 defined by 

d^(p d 2 -p 2 f (3) 
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is not easily measured precisely due to the very long transients close to the critical 
point (i.e., due to the smallness of 5) and due to finite size effects. The latter 
are absent in our simulations due to the very large lattice size. Nevertheless, 
extrapolating the data from fig.2 to t — > oo gave only a crude estimate (3 = 0.272 ± 
0.006 which is however in perfect agreement with DP where (3 = 0.2766 ± 0.0003 
31 M, M. 



In order to measure an independent critical exponent, we made in addition runs 
with initial single site damage on much smaller lattices [L < 7000) and for shorter 
times (t < 40000). Apart from the damaged sites, the initial configurations were 
randomly chosen active states (they were set to the final configuration of the first 
lattice in the preceding run by setting the z-th word to if = 0, and to -1 if 
s- 1 ^ = 1). From universality with DP we expect that at the critical point 

d~t v , t) = 0.314 ±0.001 , (4) 

which is nicely fulfilled. Off the critical point we should have 

(d) oc r 6 (f)({pi - P2)t 1/ ^ ) (full initial damage) (5) 

and 

(d) oc t v il)({jp2 — P2)t 1 ^ 11 ) (single site initial damage) (6) 

with universal scaling functions <p(z) and ip(z) which are regular at z — 0, and 
with v\\ = 1.7336 ± 0.0005. To see that our data are fully consistent with this, 
we plotted in fig.4 d/V 1 against (pf — P2)t 1 ^ for both types of initial conditions. 
We see indeed a perfect data collapse as predicted by the above ansatz. We just 
mention that similar results (again with somewhat smaller statistics and with 
significantly larger corrections to scaling) were obtained for several other values 
of pi, and allowed us to locate curve T> in fig.l with high precision. 

In conclusion we have given numerical evidence that the damage spreading tran- 
sition in the Domany-Kinzel CA is in the DP universality class, although the 
undamaged state has positive entropy. We expect this to be true in general, not 
only for the Domany-Kinzel CA. 

Of course, we have to set initial conditions such that we are not confined to 
atypical states carrying zero measure. In the present case, such atypical behavior 
would e.g. result if we would start with one of the configurations being dead (all 
Si = 0) or nearly dead (sj 7^ only in a finite region). In the latter case, we would 
then have a linear increase of d instead of eq.@. We believe that not taking 
into account this caveat is the reason why only partial universality with DP was 



observed in [g4], pq| . In these papers, 'dynamical' simulations were done where 
the active region was bounded and expanding with time. Outside this region the 
configurations were not allowed to evolve, but were (artificially) kept in atypical 
states. It seems trivial that this modification of the model can lead to violations 
of universality. 
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Unfortunately, our conjecture does not immediately apply to the case of Kauffman 
automata |33] where damage spreading had been studied quite intensively ^ . 



First of all, these models involve frozen randomness and should thus — if at all 
- be compared to DP in disordered media. Secondly, healing is not perfect in 
Kauffman models even in the phase in which damage does not spread. In this 
phase a finite damage has a non-zero probability to persist forever, and the healed 
state is not absorbing in our sense. It would be most interesting to study modified 
Kauffman models where such healing takes place (e.g. stochastic versions — apart 
from the randomness in the attribution of local rules, Kauffman models are strictly 
deterministic), and to compare them with DP in disordered media. 

We have added one more item to the already long list of possible physical realiza- 
tions of the DP transition. It is vexing that in spite of this ubiquity in models, and 
in spite of its conceptual simplicity (DP is by far the simplest critical phenomenon 
to study on a computer and to explain to a high school student), there have not 
yet been reported any experiments where the critical behavior of DP was observed 
even crudely! Maybe the present realization can lead to such an observation. 

On the more practical side, we have introduced a new and very efficient method 
for simulating damage spreading which might find applications also in similar 
problems. 

Acknowledgement: I am most indebted to Dr. M. Schreckenberg for interesting 
discussions and for pointing out the two variants of implementing damage spread- 
ing in the Domany-Kinzel CA, and to Dr. T. Lippert for help with running the 
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Figure Captions: 



Fig.l: Part of the phase diagram for the Domany-Kinzel CA. Curve C separates 
dead (left) from active (right) phases. Curve T> (which joins C at its lower end 
point, but runs otherwise entirely in the active phase) separates a healing phase 
(left) from a chaotic phase (right) where any damage has non-zero probability not 
to heal. Here the damaging is implemented according to the first variant described 
in the text. With the second variant, T> would be somewhat further to the left. The 
transition curves were determined by runs with single active/damaged initial sites, 
and demanding that the exponent r\ (see eq.(|])) has the value of DP. The precision 
of the curves is everywhere better than the thickness of the lines. Quantitatively, 



the phase diagram agrees with data from |l3|, but not with the diagram given in 



42fl . It also deviates significantly from that in fi"7 



Fig. 2: Log-log plot of the total number of damaged sites in 2016 pairs of lattices, 
with 2 22 sites each. For all curves p\ = 1, while P2 ranges from 0.3086 to 0.3155 
(from top to bottom). Initial configurations were random. 

Fig. 3: Log- log plot of the total number of damaged sites in runs where each pair 
of lattices was initially damaged at a single site. Again pi — 1. The four central 
values of pi are the same as in fig. 2. 

Fig. 4: The same data as in figs. 2 (panel a) and 3 (panel b), but plotted such that 
all data should collapse onto a single curves if eqs.(|5|||) are correct. Only data 
for t > 40 are plotted in panel a, and for t > 10 in panel b, in order to reduced 
finite-time corrections. 
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